library(tidyverse)
library(plotly)
Consider a somewhat nasty-looking curve in two dimensions. We define \(z\) to be a function of \(x\) and \(y\) that looks like this: \[
z = f(x, y) = \biggl[1 + e ^{-|x + 2y|} \cos\biggl(0.2 \sqrt{x^2 + y^2}\biggr) \biggr] \times \mathbb{I}\{-10 < x,y < 10 \}
\] In other words, it is a function that returns 0 everywhere outside of a square of width 20, centered on the origin. But within that square it returns positive values.
Here is some code that defines it:
ugly <- function(x, y) {
ifelse(x < -10 | y > 10 | x > 10 | y < -10,
0,
1 + exp(-abs(x + 2 * y) / 10) * cos(0.2 * sqrt(x^2 * y^2))
)
}
And, we can use that to visualize the surface by evaluating a lot of points and plotting it with the plotly package. Go ahead and rotate and manipulate the thing.
vals <- seq(-12, 12, by = 0.05)
pts <- expand_grid(
x = vals,
y = vals
) %>%
mutate(z = ugly(x, y))
zmat <- matrix(pts$z, ncol = length(vals), nrow = length(vals))
plot_ly(x = vals, y = vals, z = ~zmat) %>%
add_surface()
The volume under that curve is given by the definite integral: \[
\int_{-10}^{10} \int_{-10}^{10} f(x, y) dxdy = \int_{-10}^{10} \int_{-10}^{10} \biggl[1 + e ^{-|x + 2y|} \cos\biggl(0.2 \sqrt{x^2 + y^2}\biggr) \biggr] dxdy
\] That looks a little hairy. It may have an analytical solution, or it could be evaluated numerically with high accuracy. But, if you rewrite that integral as an expection over a two-dimensional uniform variable on the square from -10 to 10 in \(x\) and also from -10 to 10 in \(y\), then you get: \[
20 \times 20 \times \int_{-10}^{10} \int_{-10}^{10} f(x, y) \frac{1}{20}\frac{1}{20}dxdy =
400 \times \mathbb{E}[f(x, y)]
\] where the expectation is taken over a 2-D uniform variable from -10 to 10 in both x and y.
Given that as a preamble, your mission it to approximate that integral using Monte Carlo sampling.
Good luck! The answer is \(\approx 423.7\).
LS0tCnRpdGxlOiAiQm9udXMgZXhjZXJjaXNlOiB1c2UgTW9udGUgQ2FybG8gdG8gYXBwcm94aW1hdGUgYSBkZWZpbml0ZSBpbnRlZ3JhbCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShwbG90bHkpCmBgYAoKCkNvbnNpZGVyIGEgc29tZXdoYXQgbmFzdHktbG9va2luZyBjdXJ2ZSBpbiB0d28gZGltZW5zaW9ucy4gIFdlCmRlZmluZSAkeiQgdG8gYmUgYSBmdW5jdGlvbiBvZiAkeCQgYW5kICR5JCB0aGF0IGxvb2tzIGxpa2UgdGhpczogCiQkCnogPSBmKHgsIHkpID0gXGJpZ2dsWzEgKyBlIF57LXx4ICsgMnl8fSBcY29zXGJpZ2dsKDAuMiBcc3FydHt4XjIgKyB5XjJ9XGJpZ2dyKSBcYmlnZ3JdIFx0aW1lcyBcbWF0aGJie0l9XHstMTAgPCB4LHkgPCAxMCBcfQokJApJbiBvdGhlciB3b3JkcywgaXQgaXMgYSBmdW5jdGlvbiB0aGF0IHJldHVybnMgMCBldmVyeXdoZXJlIG91dHNpZGUgb2YgYSBzcXVhcmUKb2Ygd2lkdGggMjAsIGNlbnRlcmVkIG9uIHRoZSBvcmlnaW4uICBCdXQgd2l0aGluIHRoYXQgc3F1YXJlIGl0IHJldHVybnMKcG9zaXRpdmUgdmFsdWVzLgoKSGVyZSBpcyBzb21lIGNvZGUgdGhhdCBkZWZpbmVzIGl0OgpgYGB7cn0KdWdseSA8LSBmdW5jdGlvbih4LCB5KSB7CiAgaWZlbHNlKHggPCAtMTAgfCB5ID4gMTAgfCB4ID4gMTAgfCB5IDwgLTEwLCAKICAgICAgICAgMCwKICAgICAgICAgMSArIGV4cCgtYWJzKHggKyAyICogeSkgLyAxMCkgKiBjb3MoMC4yICogc3FydCh4XjIgKiB5XjIpKQogICkKfQpgYGAKCkFuZCwgd2UgY2FuIHVzZSB0aGF0IHRvIHZpc3VhbGl6ZSB0aGUgc3VyZmFjZSBieSBldmFsdWF0aW5nIGEgbG90IG9mIHBvaW50cyBhbmQKcGxvdHRpbmcgaXQgd2l0aCB0aGUgcGxvdGx5IHBhY2thZ2UuICBHbyBhaGVhZCBhbmQgcm90YXRlIGFuZCBtYW5pcHVsYXRlIHRoZSB0aGluZy4KYGBge3IsIGZpZy53aWR0aD04fQp2YWxzIDwtIHNlcSgtMTIsIDEyLCBieSA9IDAuMDUpCnB0cyA8LSBleHBhbmRfZ3JpZCgKICB4ID0gdmFscywKICB5ID0gdmFscwopICU+JQogIG11dGF0ZSh6ID0gdWdseSh4LCB5KSkKCnptYXQgPC0gbWF0cml4KHB0cyR6LCBuY29sID0gbGVuZ3RoKHZhbHMpLCBucm93ID0gbGVuZ3RoKHZhbHMpKQoKcGxvdF9seSh4ID0gdmFscywgeSA9IHZhbHMsIHogPSB+em1hdCkgJT4lCiAgYWRkX3N1cmZhY2UoKQpgYGAKCgpUaGUgdm9sdW1lIHVuZGVyIHRoYXQgY3VydmUgaXMgZ2l2ZW4gYnkgdGhlIGRlZmluaXRlIGludGVncmFsOgokJApcaW50X3stMTB9XnsxMH0gXGludF97LTEwfV57MTB9IGYoeCwgeSkgZHhkeSA9IFxpbnRfey0xMH1eezEwfSBcaW50X3stMTB9XnsxMH0gXGJpZ2dsWzEgKyBlIF57LXx4ICsgMnl8fSBcY29zXGJpZ2dsKDAuMiBcc3FydHt4XjIgKyB5XjJ9XGJpZ2dyKSBcYmlnZ3JdIGR4ZHkKJCQKVGhhdCBsb29rcyBhIGxpdHRsZSBoYWlyeS4gIEl0IG1heSBoYXZlIGFuIGFuYWx5dGljYWwgc29sdXRpb24sIG9yIGl0IGNvdWxkCmJlIGV2YWx1YXRlZCBudW1lcmljYWxseSB3aXRoIGhpZ2ggYWNjdXJhY3kuICBCdXQsIGlmIHlvdQpyZXdyaXRlIHRoYXQgaW50ZWdyYWwgYXMgYW4gZXhwZWN0aW9uIG92ZXIgYSB0d28tZGltZW5zaW9uYWwgdW5pZm9ybSB2YXJpYWJsZSBvbiB0aGUgc3F1YXJlIGZyb20gLTEwIHRvIDEwIGluICR4JCBhbmQgYWxzbyBmcm9tIC0xMCB0byAxMCBpbiAkeSQsIHRoZW4geW91IGdldDoKJCQKMjAgXHRpbWVzIDIwIFx0aW1lcyBcaW50X3stMTB9XnsxMH0gXGludF97LTEwfV57MTB9IGYoeCwgeSkgXGZyYWN7MX17MjB9XGZyYWN7MX17MjB9ZHhkeSA9CjQwMCBcdGltZXMgXG1hdGhiYntFfVtmKHgsIHkpXQokJAp3aGVyZSB0aGUgZXhwZWN0YXRpb24gaXMgdGFrZW4gb3ZlciBhIDItRCB1bmlmb3JtIHZhcmlhYmxlIGZyb20gLTEwIHRvIDEwIGluIGJvdGggeCBhbmQgeS4KCkdpdmVuIHRoYXQgYXMgYSBwcmVhbWJsZSwgeW91ciBtaXNzaW9uIGl0IHRvIGFwcHJveGltYXRlIHRoYXQgaW50ZWdyYWwgdXNpbmcgTW9udGUgQ2FybG8gc2FtcGxpbmcuCgpHb29kIGx1Y2shICBUaGUgYW5zd2VyIGlzICRcYXBwcm94IDQyMy43JC4K